Basic R for economics

Ignore

• SPSS
• Eviews
• MATLAB
• Stata

Try to ignore

• SAS
• Excel/VBA

Focus on

• R
• Python

Keep an eye on

• Julia

Important R people

• Hadley Wicham (tidyverse)
• Matt Dowle (data.table)
• Max Kuhn (CARET)

How to follow R

• follow rbloggers on Facebook
• follow #rstats on twitter
• stackoverflow/crossvalidated
• AalborgRUG (meetup)

compiled vs. interpeted languages

• R created in 1993 from the language S
• Quickly adopted by statisticians
• inspired pandas in python
• split into 3 dialects (base, datatable, tidyverse)

Advantages

• It is open source, and free in contrast to the expensive commercial solutions (SAS, Stata, SPSS, EViews, Matlab)
• Since it's free it's easy for others to replicate results
• A huge library of packages with tools for statistical modeling, machine learning, visualizing and manipulating data.

file extensions: .r, .rdata, .rda, .rmd

Which language?

• Python - If you want to work with engineers
• R - If you want to work with statisticians

• Vector The most basic structure, evel a single value is a vector in R. It takes different forms, the most usefull are numeric, character, factor and date.

• Matrix You can bind numeric vectors of same length together to form a matrix

• Data frame You can bind different types of vectors together to form a data frame

• List You can put all the above in a list, it can have objects of different dimensions

class(1)

[1] "numeric"

class(c(1,2,3,4,5))

[1] "numeric"

class(c("A","B","D"))

[1] "character"

class(as.matrix(c(1,2,3,4,5)))

[1] "matrix"

class(as.data.frame(c(1,2,3,4,5)))

[1] "data.frame"

class(as.list(c(1,2,3,4,5)))

[1] "list"

Very usefull in many cases

• subsetting (filtering)
• ifelse
• group_by/facet_wrap

1 > 2

[1] FALSE

1 == 1

[1] TRUE

a <- c(2,3,4) < c(1,5,6)
a

[1] FALSE  TRUE  TRUE

sum(a)

[1] 2

• Extremly important to understand
• Long (tidy data) is always prefered, but sometimes in the process you have to go back and forth between them
• In Excel this is known as pivot tables

In the tidyverse this has gone throgh a few names

• melt -> gather -> pivot_longer
• cast -> spread -> pivot_wider

relig_income[,1:4] %>% 
  head()

# A tibble: 6 x 4
  religion           `<$10k` `$10-20k` `$20-30k`
  <chr>                <dbl>     <dbl>     <dbl>
1 Agnostic                27        34        60
2 Atheist                 12        27        37
3 Buddhist                27        21        30
4 Catholic               418       617       732
5 Don’t know/refused      15        14        15
6 Evangelical Prot       575       869      1064

long <- relig_income[,1:4] %>% 
  pivot_longer(names_to = "income", 
               values_to = "value", 
               -religion)
head(long, 18)

# A tibble: 18 x 3
   religion           income  value
   <chr>              <chr>   <dbl>
 1 Agnostic           <$10k      27
 2 Agnostic           $10-20k    34
 3 Agnostic           $20-30k    60
 4 Atheist            <$10k      12
 5 Atheist            $10-20k    27
 6 Atheist            $20-30k    37
 7 Buddhist           <$10k      27
 8 Buddhist           $10-20k    21
 9 Buddhist           $20-30k    30
10 Catholic           <$10k     418
11 Catholic           $10-20k   617
12 Catholic           $20-30k   732
13 Don’t know/refused <$10k      15
14 Don’t know/refused $10-20k    14
15 Don’t know/refused $20-30k    15
16 Evangelical Prot   <$10k     575
17 Evangelical Prot   $10-20k   869
18 Evangelical Prot   $20-30k  1064

One of the most key things in R, and one of the big reasons for it's success. it's the main ingredient in pandas (python), dplyr and data.table. The operations are possible in base R, but more complicated.

We use the long dataframe of relic_income from before:

long <- relig_income %>% 
  pivot_longer(names_to = "income", 
               values_to = "value", 
               -religion)

long %>% group_by(religion) %>% summarize(n = n()) %>% head(10)

# A tibble: 10 x 2
   religion                    n
   <chr>                   <int>
 1 Agnostic                   10
 2 Atheist                    10
 3 Buddhist                   10
 4 Catholic                   10
 5 Don’t know/refused         10
 6 Evangelical Prot           10
 7 Hindu                      10
 8 Historically Black Prot    10
 9 Jehovah's Witness          10
10 Jewish                     10

long %>% group_by(religion) %>% summarize(n = sum(value)) %>% head(10)

# A tibble: 10 x 2
   religion                    n
   <chr>                   <dbl>
 1 Agnostic                  826
 2 Atheist                   515
 3 Buddhist                  411
 4 Catholic                 8054
 5 Don’t know/refused        272
 6 Evangelical Prot         9472
 7 Hindu                     257
 8 Historically Black Prot  1995
 9 Jehovah's Witness         215
10 Jewish                    682

long %>% group_by(income) %>% summarize(n = sum(value))

# A tibble: 10 x 2
   income                 n
   <chr>              <dbl>
 1 $10-20k             2781
 2 $100-150k           3197
 3 $20-30k             3357
 4 $30-40k             3302
 5 $40-50k             3085
 6 $50-75k             5185
 7 $75-100k            3990
 8 <$10k               1930
 9 >150k               2608
10 Don't know/refused  6121

y <- rbinom(100, 1, prob = 0.5)
x <- abs(rnorm(100)) + y

fit <- glm(y ~ x, family = "binomial")
fit$coefficients

(Intercept)           x 
  -4.805890    4.082138 

ggplot(tibble(x,y), aes(x=x, y=y)) +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE, color="grey") +
  geom_point() + labs(title="Logistic Regression")

For-loops are easy and usefull, but can sometimes be slow. I calculate simulate data and calculate a logistic regression, i then run this 10 times in a for-loop each time saving the coeffecient.

fits <- NA

for (i in 1:10) {

  y   <- rbinom(100, 1, prob = 0.5)
  x   <- abs(rnorm(100)) + y * 2
  fit <- glm(y ~ x, family = "binomial")

  fits[i] <- fit$coefficients[[2]]

}

kable(fits)

• operators: <, <=, ==, >=, >, %in%, !
• table()
• subset()
• ifelse()
• simulate: rnorm(), runif(), rbinom()
• class(), str(), summary()
• seq(), rep()

df <- data.frame(x = sample(c("A","B"), 
                            prob = c(0.2,0.8), 
                            size = 100, 
                            replace=T),
                 y = rep(c("Apple","pear",
                           "Apple","bananna"), 
                         25))
head(df)

  x       y
1 B   Apple
2 B    pear
3 B   Apple
4 B bananna
5 B   Apple
6 B    pear

table(df$x, df$y)

    Apple bananna pear
  A     4       4    4
  B    46      21   21

Lets say we have a more normal mathematical function, we can write it in R

\[ f(x) = \bigg(x^3cos\frac{x}{2}+\frac{1}{2}\bigg)\sqrt{4-x^2} \]

fun <- function(x) {
  (x^3 * cos(x/2) + 1/2) * sqrt(4 - x^2)
  }

fun(1)

[1] 2.386043

fun(-1)

[1] -0.6539922

What if we would like now to integrate this function within an interval?

\[ \int_{-2}^{2}\bigg(x^3cos\frac{x}{2}+\frac{1}{2}\bigg)\sqrt{4-x^2}\ dx \]

integrate(fun, -2, 2)

3.141593 with absolute error < 2e-09

ggplot(data.frame(x = 0), aes(x)) + 
  stat_function(fun = fun) + 
  xlim(-2,2)

Question A (hard)

With 1000 coin tosses you get a head 530 times \( p=0.53 \), is this a fair coin? calculate the confidence 95% confidence intervals with the following formula:

\[ CI=\pm 1.96 \sqrt{\frac{p(1-p)}{n}} \]

How can you know if this formula is true? We can calculate the coverage and verify that the calculation of confidence intervals.

1. Simulate 1000 coin tosses with \( P = 0.5 \) and calulate UCL and LCI (use rbinom() to simulate cointosses)
2. We know that 0.5 should be within this interval 95% of times, so with logical operators check wether 0.5 is contained within the confidence intervals
3. In a forloop run the above 100 times and each time save the TRUE or FALSE value within a vector
4. sum the vector

Question B (easy)

You simulate a x and y vector each with 1000 obs in R with rnorm(), and calculate

set.seed(1)
x <- rnorm(1000)
y <- rnorm(1000)
z <- rnorm(1000)
lm(y ~ x + z)

Call:
lm(formula = y ~ x + z)

Coefficients:
(Intercept)            x            z  
  -0.016544     0.005328     0.022451  

check the calculation of the coeffecients using the matrix formula

\[ \hat{\beta}_{OLS}=(X'X)^{-1}X'Y \]

If the results don't fit, did you forget intercept? lm(y~0+x)

Inequality simulation

Imagine you are offered the following bet with the probarbility of 50% \( P \) to either increase or decrease your total wealth \( M \)

\[ result = P(1.5*M) | P(0.6*M) \]

1. What is the expected value of this bet \( E(result) \)

2. Since this is a good bet, what would the effect be on the wealth of soeicity if we let all 5 million people in denmark do this bet 100 times? create

3. What would be the effect on the average person?

See code on the right for help:

set.seed(3)
individuals <- 100
bets        <- 100
cash        <- 100

mat <- matrix(data = NA, 
              nrow = individuals, 
              ncol = bets)

for (j in 1:individuals) {
  for (i in 1:(bets-1)) {
    cash[i+1] <- ifelse(rnorm(1) > 0, 
                        cash[i] * 0.6, 
                        cash[i] * 1.5)
  }
  mat[j,] <- cash
}

m1 <- lm(wage ~ education + 
                experience + 
                age + 
                gender, 
         data = CPS1985)

m2 <- lm(wage ~ education + 
                experience + 
                age + 
                gender + 
                ethnicity, 
         data = CPS1985)

m3 <- lm(wage ~ education + 
                experience + 
                age + 
                gender + 
                ethnicity + 
                married, 
         data = CPS1985)

stargazer(m1, m2, m3, 
          column.labels = c("Good", 
                            "Better", 
                            "Best"),
          single.row = TRUE)

CPS1985 %>% 
  ggplot(aes(age, wage)) + 
  geom_point() + 
  geom_smooth(method = "lm") +
  labs(title="Does age influence wage?")

CPS1985 %>% 
  ggplot(aes(age, wage)) + 
  geom_point() + 
  geom_smooth(method = "lm") +
  facet_wrap(~ occupation, nrow = 2)

library(sf)

geo_sf <- read_sf("https://dawa.aws.dk/Sogne/?format=geojson")
KMGALDER <- sdk_retrieve_data("KMGALDER", Tid = "2019", Køn = "TOT")
KM <- subset(KMGALDER, KMGALDER$KØN != "Total")
KM$kode  <- str_extract(KM$SOGN, "[:digit:]+")
KM$INDHOLD <- KM$INDHOLD/10

KM %>% 
  left_join(geo_sf, ., by = c('kode' = 'kode')) %>% 
  ggplot() + 
  geom_sf(aes(fill = INDHOLD, color = INDHOLD)) +  
  coord_sf(xlim = c(8.5, 11.5), ylim = c(56.5, 57.5), expand = FALSE) +
  theme_minimal()

• Create user to get API key

library(Quandl)
library(tidyverse)
library(lubridate)

ZERO <- Quandl("FED/SVENY", api_key = key)
head(ZERO)[1:4]

        Date SVENY01 SVENY02 SVENY03
1 2019-11-29  1.6760  1.6416  1.6291
2 2019-11-27  1.6784  1.6520  1.6405
3 2019-11-26  1.6448  1.6118  1.5971
4 2019-11-25  1.6506  1.6282  1.6199
5 2019-11-22  1.6508  1.6296  1.6223
6 2019-11-21  1.6364  1.6169  1.6114

range(ZERO$Date)

[1] "1961-06-14" "2019-11-29"

difftime(range(ZERO$Date)[2],range(ZERO$Date)[1])

Time difference of 21352 days

nrow(ZERO)

[1] 14577

Quandl("FRED/BASE", api_key = key) %>% 
  ggplot(aes(Date, Value)) + 
  geom_line() + 
  geom_ribbon(aes(ymin=0, ymax = Value), alpha = 0.4) +
  labs(title = "US high-powered money (monetary base)",
       subtitle = "What happened in 2008?\nWhy didn't this result in inflation?")

Quandl("FRED/WRBWFRBL", api_key = key) %>% 
  ggplot(aes(Date, Value/1000)) + 
  geom_line() + 
  geom_ribbon(aes(ymin=0, ymax = Value/1000), alpha = 0.4) +
  labs(title = "Liabilities and Capital: Other Factors Draining Reserve Balances:",
       subtitle="Reserve Balances With Federal Reserve Banks: Wednesday Level")

We estimate a linear model with an endogeneity problem

\[ y_i = 0.1 W_i + \varepsilon_i \]

\[ W_i = e^{-x^2_i}+u_i \]

\[ \varepsilon_i,e_i ~ N (0,\Sigma) \qquad \Sigma = \begin{pmatrix} 1 & 0.9 \\ 0.9 & 1 \end{pmatrix} \]

library(mvtnorm) # Multivariate rnorm
library(gmm)     # Generalised Method of Moments

n <- 1000
e <- rmvnorm(n, sigma = matrix(c(1.0, 0.9, 0.9, 1.0),2, 2))
cor(e)

          [,1]      [,2]
[1,] 1.0000000 0.8957488
[2,] 0.8957488 1.0000000

x <- rnorm(n)
w <- exp(-x^2) + e[,1]
y <- 0.1 * w + e[,2]
df <- tibble(y, w)

ggplot(df, aes(w, y)) + 
  geom_point(shape = 21, alpha=0.2) + 
  labs(title = "Data with an endogeneity problem")

CARET

is an standardized implementation of algorithms for prediction, these include:

• Regularized regressions (Elastic Net) glmnet
• Linear Discriminant Analisis (LDA) lda
• Support Vector Machines (SVM) svmLinear
• Regression Trees (bagging/boosting) tree
• Random Forrest rpart randomforrest rf
• Neural Nets (Keras for more advanced) nnet

Max Kuhn that is behind CARET is working on a new and improved framework tidymodels (parsnip) but this is still a work in progress.

library(caret)

index <- createDataPartition(CPS1985$gender, p = .75, list = FALSE)

training <- CPS1985[index,]
testing  <- CPS1985[-index,]

fit.net <- train(gender ~ ., 
                 data   = training,
                 method ='glmnet')

pred <- predict(fit.nnet, newdata = testing)
table(pred, testing$gender)
#confusionMatrix(data = pred, reference = testing$gender)

Key bindings